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APS TRACT 

There are many I'r-actica] sltuatjons where twu parts are 
direct iy in contact subjected tc* dynamic or varying lt>ad. TTiis 
leads to development of high stresses at contact region and may 
become a cause of failure. Due to variation in the load, the 
contact area and contact conditions keep changing and are not 
known a priori at any load. The problems 'are non-linear. The 
regions of slip and no-slip occur at tlie contact .surface. The 
conditions of contact are determined by t,lie kinematic con.straints 
and Coulomb’s law of friction. Finite element analysis of ^ two 
dimensional elastic quasi-static and roiling contact problems is 
presented. As tiie load varies with time, direct application of 
finite element analysis i.s difficult as it leads to non-unique 
solutions If load steps are large and mesh Is coarse. In the 
present analysi.s, tlie principle of minimum dissipation of energy 
is used along with finite element analysis. 

The increoiental loading and unloading is emplc:'yed in 
quasi -static contact problem and tlie concept of travelling finite 
element.s is made use of in solving tiie steady rolling contact 
problem. The combined incremental and iterative procedure is 
adopted to solve the.se problems. The displacement function based 
finite element method is used to Implement above algoritlira. Tlie 
results obtained aro in good agreement with physical reasoning. 
The computational time and efforts have been greatly reduced by 
using skyline-storage for global stiffness matrix and 
sub-. structuring metliod.s . 



CHAPTER I 


INTRODUCTION 


l.J GENERAL INTRODUCTICHW 

There are many instances in engineering structures, wherein 
the components are contacting one another and transfer external 
forces through the contact surfaces, such as shrunk-fit shaft ai^d 
rotor, gear teeth in mesh, flanged connections of pressure 
vessels, dovetails of turbine blades, rolling wheel, etc. Hig}i 
stresses occur near contact zone leadingtoci'ack initiation and 
failure. Tliis situation may occur in almost all types of 
mechanical assemblages. Therefore, it is important for the 
designer to understand the stress and deformation pattern in the 
vicinity of the contact surfaces. Accurate stress and strength 
analysis of such structures requires careful understanding of the 
frictional contact conditions. 

In practical situations, the contact problems of different 
natures are encountered. The present investigation is aimed at 
elastic bodies in contact. These problems can broadly be 
classified into three categories: 

A: Stationary Contact with no Friction 


This occurs when two smooth bodies come into contact with 
each other along plane surface and transfer forces through it. No 
friction develops here. This is a linear problem because the 
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contact area remains unchanged during loading and also reversible 
as the system is conservative. This is the most simple problem 
and can be solved by direct application of principles of 
elasticity. 


B. Varying Contact With no Friction 


When two curved bodies with smooth surfaces come into 
contact, the area of contact is a function of load applied. Wlien 
the contact area increases with the load, it is termed as 
advancing contact and receding contact vice versa. The change in 
contact area with loads renders the problem non-linear. However, 
it is still of reversible nature due to absence of 
non-conservative forces. Iterative procedures are adopted for 
solving such problems. 


C- Varying Contact with Friction 


When the nature of surfaces are frictional and the contact 
area varies with load, the problem comes under the category of 
varying contact and non-linear. The presence of friction mali.es 
the system irreversible and thus making it complex. This, in 
fact, is the most general case of contacting problems. 

The present work is an attempt to solve these types of 
problems. Here, under the action of external loads, the domain of 
contact area and distribution of tractions are unknown over the 
contact region, thus, the contact problem Is a non-lineai' problem 
with unknown boundary conditions. The existence of friction may 
lead to the formation of regions of slip and no-slip within the 
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contac-b surface. The relative Diagnitudes of tangential and normal 
tractions governs the status of contact surface. Also, the energy 
dissipated by the frictional forces makes the problem 
irreversible. Due to this the loading history also influences the 
final stress pattern. 


Many analytical and numerical approaches have been attempted 
by many investigators in solving contact problems. The present 
work is an attempt towards using finite elenient method and 
principle of minimum dissipation of energy in solving two 
problems. 


The first one, deals with a large plate with a circular hole 
having an inclusion subjected to slowly vai-ying cyclic loading. 
This problem is of practical interest as all machines con.sists of 
many plates and are subjected to vibration. Another kind of 
commonly encountered contact problem is a rolling contact problem, 
where a body, having curved surfaces or boundary, rolls over other 
body which may be plane or curved. Such a phenomenon is found in 
many practical situations like tire-road contact, rail-wheel 
contact, ball and roller bearings, etc. In all these cases, 
during the rolling motion of one body over the other, the loads 
are transferred through contacting surfaces which keep changing 
continuously, 

1-2 LITERATURE SURVEY 

Contact problems have long been of considerable interest. 
Since, Herts (1882) gave his well known classical solution of 
frictionless contact problem, further developments have taken 
place. The need for this in-depth study has come from practical 
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considerations particularly with respect to influence of 
interfacial friction. The following three basic approaches were 
adopted in solving these problems. 

The earliest solutions have been obtained using integral 
equation methods. Existing exact solutions of contact stresses 

are the product of sophisticated mathematical analysis for 
idealised configurations. In many real situations, it is not 
possible to find suitable models for which closed form solutions 
are available. Tlie exact solutions have been contributed by 
Mukhelishvili [1963] and Gladwell [1973], etc. 


In later development, contact problems were considered as a 
special case of constrained optimization of either total potential 
or complementary potential energy. The minimisation was 
formulated as a mathematical programming problems [ Hung, 1980 ] . 

In the third approach, contact conditions were imposed 
directly from kinematic considerations by imposing geometric 
compatibility of contacting surfaces during the incremental 
loading process. The main advantage of this metfiod is that finite 
element procedures can be effectively used and the friction 
conditions can be imposed with ease. Finite element formulation 
can also be categorised in three types. 

The flexibility method was developed by Francavilla and 
Zienkiewicz [1975] for frictionless elastic contact problems and 
then extended to frictional contact by Sachdeva and Eamakrishnan 
[1981]. Jing and Liao [1990] introduced the concept of relative 
displacements and contact double forces. Total loading and 
compatibility of displacements was used to determine the contact 
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forces. In second category, a stiffness aatrix was used in the 
finite element formulation. Wilson and Parsons [19703 used 
differential displacements and a load incremental process to treat 
the frictionless coxitact problems. Chan and Tuba [1971] and Ohte 
[1973] studied the frictional contact problems for plane and 
axisymmetric elastic problems. Fredrichsson [1976] proposed a 
constitutive model to account for friction at contact surface. 
Gaertner [1977] used an experimental friction law and a modified 
element to treat the connecting rod and Herts problem. Okaisioto 
and Nakajawa [1979] gave a finite element solution using 
incremental load theory. Torstenfelt [1983] and Rahman, Cook and 
Rowlands [1984] used Iterative procedure for analysis of 
elastostatic contact problem with friction. Bathe and Choudliary 
[1985] solved for frictional contact under large deformation. 
Other researchers who contributed in the elastostatic field are 
Chandrasekaran, Haisler and Goforth [1987]. A variational 
approach, where the e<iuilibrium and the boundary conditions were 
put into variational form, was adopted by Oden and Fires [1983]. 
They used a non-classical, non-local friction law and variational 
principle under which these laws hold. 

Importance of rolling contact has also been considered 
by many authors. The different analytical and numerical solutions 
have been given by Poritsky [ 1950 ], Johnson [ 1958 ], Kalker [ 
1979 ], etc. Kalker [1979] used the variation principle of Duvout 
and Lions for dry friction in solving three dimensional rolling 
contact problem with Hertzian normal contact. Liu and Paul [1989] 
developed a program to solve three dimensional rolling contact 
problems for arbitrary contact patches under arbitrary pressure 
distributions. The problem was considered with friction and 
non -Hertzian Pressure distribution. Fadovan, Tovichakchaikul and 
Zeid [1981,1984] developed travelling finite elements based on the 
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use of a moving total /updated Lagrangian observer. They included 
large deformation kinematics and kinetics in the formulation. 


The rolling contact problem was also described by a 
variational inequality (hyperbolic/elliptic) and solved 
numerically using incremental finite element method by Zochowski 
and Mysliriki [1991]. They proposed a quasi-static approach and 
regularization of friction conditions. The problem was formulated 
as optimization problem and solved for normal contact tractions 
initially. The tangential contact stresses were calculated from 
regularised friction conditions. 


The present work on the rolling contact, is an attempt to 
solve the problem by finite element method without assuming any 
distribution of tractions on the contact patch. The friction 
conditions are also not regularised and law of Coulomb’s friction 
is followed. The concept of travelling finite elements is used to 
simulate the rolling effect. 


Either in quasi-static or steady rolling problems, the 
friction contact problem has multiple solutions satisfying 
equilibrium conditions. As suggested by Kishore.S.H. [1979] and 
mentioned in the context of plasticity problems [ Jolmson and 
Mellor, 1973 3, it is proposed to use the concept of minimum 
dissipation of energy. That is, out of all the possible 
equilibrium configurations satisfying the new load conditions, the 
actual one leads to minimum dissipation of energy from the earlier 
state. 
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1.3 SCOPE W PRESENT WORK 

The author studied the steady rolling contact ai^d 
quasi-static contact problems. In present study, finite element 
method coupled with principle of minimum dissipation of energy has 
been employed. In the rolling contact problem, the concept of 
travelling finite elements is used and a simple algorithm has been 
presented. Ko assumption regarding the distribution of tractions 
on the contact patch is made. In quasi-static contact problem, an 
incremental loading process has been adopted. Tlie continuous 
loading and unloading, ranging from given maximum tension to that 
in compression, has been applied and corresponding response is 
obtained. Computer programs have been developed to solve these 
two types of contact problems and example problems are solved. In 
quasi-static contact problem, a large elastic plate with a hole 
having an oversized rigid inclusion is investigated. In rolling 
contact, a problem of steady rolling of a cylinder on a rigid 
plane surface, is solved. 

The description of the problem formulation and solution 
techniques are detailed out in the following chapters. The theory 
of contact problem and the principle of minimum dissipation of 
energy are discussed in chapter II. Chapter III deals with the 
development of mathematical model and the solution scheme 
employed. The discussion on results is in chapter IV and chapter 
V summarises conclusion of the present work and suggestions for 
future investigation. 



CHAPTER II 


THEORY OF CONTACT PROBLEMS AND 
MINIMUM DISSIPATION OF ENERGY APPROACH 


2.1 BASIC THECtfiy OF CC»4TACT PROBLEM 

Consider two generic bodies which are arbitrarily denoted as 
contactor and target ( fig. 2.1a ). The bodies are brought into 
contact by prescribed tractions and displacements applied on the 
surface S+ and S respectively. The target body caii be either 
rigid or deformable. On coming into contact, the contact forces 
develop in the region of contact ( fig. 2.1c ), where neither the 
contact pressure nor the area of contact is known a priori. 

The basic condition,s of contact along the contact surface is 
that no material overlap can occur ( fig. 2. lb ) and as a result 
contact forces are developed that act upon the target and 
contactor. The normal tractions can be only compressive in nature 
and the tangential tractions satisfy the law of frictional 
resistance. According to Coulomb* s law of friction, the relative 
motion between two adjacent particles on the contactor and the 
target in contact can not occur, as long as jrj < p > where 

is the compressive traction, t is the tangential traction and p 
is a coefficient of friction. Wlien Jt} > the particles 

slide over each other and the motion continues as long as the 
tangential traction is developed to equal the limiting friction 
force Thus, the regions of no-slip and slip are identified 

and the following conditions are imposed as given below ”• 


Ci3 stick or no-slip conditions 
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F® 
n 

( 2.1 ) 

Jk S 

= Normal displacements of adjacent particles of body A 
and B, which are in contact. 

j = Tangential displacements of adjacent particles of 

body A and B, which are in contact. 

F.^^, F^ are corresponding normal and tangential contact forces 


A 

F^ =- 
n 

< =- 


C i v> si ip c ondi -t i ons 


u 


n 




F* =-F® 
n n 

F* =-F® 

^t ^t 

The direction of 
opposite in nature. 


( 2.2 ) 

and the tangeiitial displacement are 


The same conditions of contact and friction apply to both 
quasi-static and rolling contact problems. 

Wlien a roller moves on a rigid plane surface without sliding, 
the contact develops between roller and ground. Within this 
contact area, some regions may be in slip and the other in no-sllp 
condition. Generally, the micro-slip region is found in the 
leading edge and macro-slip occurs over trailing edge [Moore, 
1975]. Since slip Is an irreversible process, different loading 
histories would lead to different patterns of deformations- Three 
cases of rolling conditions are found in practical situations. 
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viz., free, driven and braked. In free rolling i.e. rolling under 
the action of a purely normal load, the extent of slip is small 
whereas when the contact surfaces transmit an externally applied 
tangential force or twisting couple, an appreciable amount of slip 
occurs . 

Uniform rolling motion is a steady process in which the loads 
and the fields of stress, strain and deformation do not change 
with time with respect to an observer translating with the speed 
of roller. This is in contrast to quasi-static loading, where the 
stress-strain field and deformations changes with load. Therefore, 
the rolling contact problems are also termed as moving contact 
field or travelling load problems. 


a.E MINIMUM DISSIPATION OF ENERGY APPROACH 

The elastic contact problem with finite friction is an 
irreversible prc>blem, the solution to which depends upon the 
loading history. Here, not only the value of load but also the 
history of loading decides the direction of frictional resistance 
and hence the slip direction. In this problem, the loading ond 
unloading are done gradually and the process is assumed to be 
quasi -static. 

To explain the minimum dissipation approach, consider a 
system shown in the fig ( 2.2a ). Let P the load applied and F 
be the friction force developed. The force balance relation,ii 

K X + F = P ( 2.3 ) 


and 

a- Lso F is constrained as follows 


F s lujNl 


( 2.4 ) 
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The above relation and constraint has a feasible range for 
the displacement X for load level PI and P2, shown in fig. ( 2.2b 
), is given by broken line. If the system were at a, at the load 
level PI, due to the load Increment, it is proposed that it would 
go to the point b wherein, dissipated energy is minimum, in this 
case sero, and is feasible point. If it were at the point a’ , at 
the load level PI, even then it would move to the point b since it 
is that feasible point which has minimuni dissipation, but unlike 
the previous case there is a finite amount of dissipation in this 
case. 


The approach adopted by the present work is that after 
establishing the total contact area, the target is to decide about 
slip and no-slip regions. The principle of minimum dissipation 
of energy is employed to determine the extent of tangential 
displacement in the slip region. 

The procedure of minimi aati on of energy of dissipation and 
the search for the feasible state is formulated as a constrained 
minimization problem, with the 

dissipated energy expression, as the objective function 

together with a penalty due to constraints, e. The expression for 
is given as : 


m 

*d = 2 1 ‘"k 

k=l 


:<c ((F^ + (F^ )^) / 2 I 

k k 


( 2.5 ) 


where the summation is over all m-contact node, denote the 
tangential displacement of the node k on the contact surface at 
i"^^ load step and wl .. 


n. 


represent the tangential and normal 


. th 


nodal force at i load step. 

The penalty term « , takes 

conditions given by : 


into account the friction 
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2 

+ C ^ j (ct^ 


* CjF^j - 1F^1>^ 

( 2.6 ) 


The first term in above expression is considered only for 
those nodes where j F^l 2: #j) F^j . X & C are penalty parameter and 
weighing factor, respectively . 

To obtain the solution of constrained minimisation problem, 
the objective function is modified by introducing penalty terms 
for the constraints. Now, the modified objective function is 
minimized, using a suitable unconstrained optimization technique. 

Thus, the objective function to be minimized is '. 

+ « ( 2.7 ) 

When ’*> is minimized, « usually becomes zero showing that all 
the nodes satisfy contact conditions and also assumes a minimum 
of such possible states, then that state is characterized as a 
feasible state. 

With the above notation, the sticking status of any node is 
represented as follows: 

When ^ 

e = 0 and ® 

When e is not equal to aero then the set of 
unknowns are found by minimizing Thus, the extent of slip 

is determined. 
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Prescribed 



displace men-ts 

(a) Condition prior to contact 



(b) Condition at contact 


^forces 



Applied external 
forces 


Contact 


Contact 

surfaces 



(c) Forces acting on contactor and target body 


Fig. 2.1 Schematic representation of a contact 
problem . 




Fig.2.2(a) Example contact configuration. 



Fig.2.2<b) Feasible region. 



CHAPTER III 


FINITE ELEMENT METHOD AND SOLUTION ALGORITHM 


The solution of the contact probleia is determined in two 
stages. The first part is to determine the area of contact and 
the second paz't is to find the feasible state by principle of 
minimum dissipation of energy. The second part is the computer 
time intensive and hence sub-structuring methods are used to 
simplify the repeated calculations. 

3.1 FINITE ELEMENT FORMULATION 

To obtain a finite element solution, the following 
assumptions are made : 

1. The materials are linear elastic. 

2. The deformations are small. 

3. Plane strain condition prevails. 

4. The body forces are negligible. 

5- Tlze frictional force acting at the contact surface follows 
the Coulomb' s friction law. 

6. In case of rolling contact problem, the rolling speed is 
constant . 

The rolling speed and mass are small enough, so that inertia 
forces can be neglected. Hence, the problem is quasi-static. 


7 . 
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Consider a two-dimensional elastic contact problem with 
friction and separation. It is decided to use the stiffness 
formulation using displacement functions^ as the contact algoritlim 
can be implemented effectively and frictional condition can be 
enforced with ease. 


The total potential energy functional in terms of stress 
functions is given as : 

" = I / < «ij > ‘i'' - / “l> '3'^ - j ' ■^1 is 

V V s, 

{ 3.1 ) 

where , 

c’'. . , €. . stress and strain tensors respectively 

-X- tj 


Xi.Ui 


Qi^Ui 


body force components and corresponding 
displacements in volume domain V. 

surface tractions and corresponding 
displacements over the surface traction boundary 


The first term in the above expression is strain energy and 
the last two are the work of body forces and surface tractions 
respectively. 


To derive finite. element formulation for the given 
functional, the domain is discretized into standard finite 
elements. Since the highest order of the derivatives present in 
the functional is one, minimun:i C° continuity is essential for the 
shape functions in the element selected. The simplest element 
satisfying this requirement is a constant strain triangular (CST) 
element i.e. a 3-noded triangular element. The derivation of 
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element stiffness matrix is described briefly in the following 
paragraph : 


Assume displacement functions for 3-noded triangular element 

as : 


u = + agX + a^y 

u^ = + b2X + b^y 

( 3.2 ) 

where u^ &. are displacement field in x &. y directions, 
a' s and b' s are unknown coefficient. 


The displacement functions can be represented in teinns of 
nodal degrees of freedom as '• 


u^ = N . + N . u^ + N.> u^ 

11 J J k k 




( 3.2 ) 


where N' s are interpolation function for CST elements, given 


as : 


where 




1 

2A 





X + c 





( 3.3 ) 
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= ^k ■ 


similarly, and Nj^ can be determined. 


The strain components for the element 
displacements are given as : 


in terms 


{ s } ^ 


du /^x 
&vr/»y 


du^/<?y + du^/3x 


{/.} = [ B 3 { u } 


where , 


[B] = 


( 3.4 ) 


{u} = 
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The stress-strain components are related by the material 
property vector [D] as; 


{a} = [D] {e} 


( 3.5 ) 


For Plane Stress, 


CD] 


For Plane Strain, 


E 


1-v' 


0 


CD] = 


E 


1 

0 


0 

0 


(l+j>)(l-2v) 


1-jj 

l-i> 

0 0 


0 

0 

l-2i> 


( 3.6 ) 


( 3.7 ) 


where E and v are the modulus of elasticity and Poisson' s 
ratio. 


The given functional, neglecting body forces : 


" I / ‘ '■ij *13 > 'i'' - X < H “i > ■'S 

Taking the matrix form and by substitution - 
n = 1 J {€}'^[D] {€} dV - J { u. ) dS 


( 3.8 ) 
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^ I J {u}'^[B3'^[D] [B] {u} dV - {u}’^{F} 

V 

= I / {u}'^CK3 {u} - {u}'^ {F} 

V 

where CK3® = / CB3'^ [D] [B3 dV 

V 


= Stiffness matrix of a triangular element. 


Minimizing n with respect to nodal displacements and it 
results into set of equations 

[K3 {u} = {F} ( 3.9 ) 

where 

[K3 = Global stiffness matrix = ^ [ K 3^ 

{u} = Nodal displacement vector 
{F} = Nodal force vector = ^ { F 


The system of equations for two bodies in contact can be 
combinedly written as 


: 0 


' 0^ ' 


F^ 

_ 0 i K® _ 

*1 

D® 

k J 

>■ z: -i 

F^ 

k 


where , 

, D® 
F^, F^ 


stiffness matrix of two bodies A & B 
respectively . 

displacement degrees of freedom of bodies A & 
B respectively. 

the force degrees of freedom of bodies A and B 
respectively 


) 
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For any problem, bhe accuracy of "the solution also depends 
upon the mesh size but at the cost of computational efforts. A 
good compromise is to consider a fine mesh in the region of higVi 
stresses i.e. in the vicinity of contact surface and coarser 
mesh away from this region. 

In the present investigation, for simplicity, it is assimed 
that one of the two contacting body is rigid, whereas the other 
one is deformable which needs to be discretized only. SyDmietry 
has been made use of wherever present. At the contact surface, 
local coordinate system, having tangential and normal directions, 
are considered. For solving above equations the contact 
conditions have to be enforced and they are described in the 
following section. 


3.2 CONTACT CONDITIONS 

As stated above, one of the two bodies in contact is taken as 
rigid one, the set of equations corz'esponding to that body need 
not be considered. Hence, eq. ( 3.10 ) simplifies to 

[ K ] { D } = { F } ( 3.11 ) 


The solution of the above set of equations is subjected to 
boundary conditions. There are two kinds of boundary conditions 
in contact problems, viz. prescribed and contact. The 
prescribed boundary conditions are in terms of externally 
applied tractions and displacements. Tlie contact boundary 
conditions are not known a priori and are determined in two 
phases . 
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Phase I ; Determination of contact area 

Initially, a set of possible nodes is considered in contact 
and the response is obtained at given load level. If at a node, 
the normal contact force turns out to be tensile then, that node 
is set free. Wiereas, if the displacement of any free node, in 
normal direction of contact surface, appears to penetrate 
geometrically or touches the rigid surface, it is included in the 
contact region. 

Thus, the following conditions are in^posed when the node 
changes their state : 

(i) Change from Contact to Separation 

In the separated condition, the node is treated as a free 
node and their nodal force components are equated to zero. 

^n = ® 

^'t = ^ 

F ,F. are normal and tangential contact forces respectively. 

Tl X 

(ii) Change from Separation to Contact 

The displacement conditions in normal and tangential 
directions are modified so that the given node becomes a contact 
node as ; 
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where , 

is a normal nodal displacement, 

^ is a clearance between the rigid surface and the 
corresponding node in normal direction before 
application of the given load 
u^ is the tangential nodal displacement which is 

negligibly small. 

Using the above Conditions at the corresponding nodes, the 
response in terms of reaction contact nodal forces at the contact 
zone and displacements of free nodes is obtained. Again, the 
conditions of contact are checked and the above procedure is 
repeated, till the contact conditions no more changes at the given 
load in successive iterations. Now, all the nodes in contact will 
have compressive normal force. 


Phase II : Determination of Slip and No-slip Zones in the 
Contact Region 

Having determined the contact area at a load level, the next 
phase is to check the friction condition to establish the regions 
of slip and no-slip as given below: 

Let ju be the coefficient of friction and and F. be the 

n t 

normal and tangential nodal contact forces. 

(i) when jF^j :S }F^j , both the contact surfaces adhere to each 
other i.e. there is no relative tangential displacement 
between the node and the rigid- surface . 
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For this case, the following conditions are imposed as: 
^n -- 

Ut = 0 

(ii) When jF^l > > "the node tends to slip over the rigid 

surface i.e. the node will have relative tangential 

displacement The nodes are allowed to slip hy extent such 
that jF^j = jF^j . The direction of F^ is chosen such that 
it opposes the tangential movernent of the node with respect 
to rigid surface. The new solution satisfying all conditions 
is determined by the approach discussed in sec. 2.2. The two 
types of problems are studied in the present work. 

3.3 APPLICATION TO QUASI -STATIC LOADING 

In this, a plate with a hole with an oversized rigid 
inclusion is studied. The contact between the plate and inclusion 
is due to mechanical friction only. The plate is subjected to 
farfield uniform loading. The loading and unloading processes are 
done gradually and after each load increment, the solution is 
updated. In such problems, the load increments are chosen in such 
a manner that at a time either one node comes into contact or just 
one node separates, and the changed contact conditions at the node 
are assumed to become operative for subsequent load increments. 
The expressions for interpolating the incremental load is 
discussed below as : 

First, solve for a nodal disp'lacement increment scad a 

contact force increment vector {AF > using the contact conditions 

n 

of previous step and a test load {AP}. {U„} and {F^} are 

n. n 

solution at load P^. 

n 
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The new load is found as 




( S.12 ) 


where a is a load increiftent.al ratio. The value of a is 
smallest of ( It is a value of load incremental ratio at 

different nodes ) calculated for the nodes violating their 
existing contact condition.^. The following two cases are 
considered : - 

1. From Separation to Contact Condition 


ivT .1-^x7 4a* 

N,n+1 N,n n 
= 0 


Al 


N,n 


( 3.13 ) 


So, 


ot* = - I „ / A I „ 

n N,n N.n 


( 3.14 ) 


and A are the clearance between two bodies in normal 
direction before the application of incremental load and change 
in clearance over the incremental step respectively. 


2. From Contact to Separation Condition 


So, 


Fx7 ^1 = Fm + a * (AF„ ) 
N,n4l N,n n N.n-^ 

= 0 


a * = - F„ / AF„ 
n N,n N,n 


( 3.15 ) 


( 3.16 ) 


F„ is normal contact force component before application of 
N 

incremental load. 
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is incremental normal contact force component 

The same criterion of contact and friction apply as discussed 
in the sec. 3.2. In this problem, because of the interpolation of 
load increments, a new node is included in the contact aone when 
it has just touched the rigid surface, hence the normal clearance 
^ is always aero. Wlien the coi'idition of slip is encountered, the 
extent of slip is determined by the approach discussed in section 
2-2, The dissipation of energy, (Eq.2.3) and penalty term, « 
(Eq.2.4) are reproduced below : 


$ 


d 


m 

y j - «^) * {(F + (F^ )^) * / 2 I 

k=l ^ ^ 

( 3-17 ) 


whei’e represents the sum of the dissipation of energy at 
all m-contact nodes. 

The penalty term « , takes into account the friction 

conditions as given by : 


m 


- 1 


fd*',. I - p 


k=l 


k 




+ C j (« 


i-1 




* CIF.I 


Thus, the objective function is : 

't = + « 

d 


HFJ>| 

( 3.18 ) 
( 3.19 ) 


The Objective function which includes the penalty, « and 
the dissipation of energy, 5^, is minimiaed by a suitable 
unconstrained minimiaatlon technique. Hookes and Jeeves pattern 
search method [Rao S.S., 1990 and Kuester-Mize, 1973 ], which is a 
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unimodal, sequential search type, is adopted in determining the 
tangential displacements, represented as the geometric variables 
such that the dissipation is minimum and all other constraints 
are satisfied. For the kind of objective fui^ction, formed in this 
problem, the Hooke and Jeeves pattern search method is found to be 
quite effective. This method can handle no. of variables 
simultaneously and convergence is guaranteed. 


In this method probing fui^ction evaluations are made Itx all 
the directions and as shown in fig. 3.2, for two variable seai'ch. 

The variables are perturbed sequentially and pattern direction, 
S, is found. The variables are given one step in the pattern 
direction before the search for a new direction is made. If 
search fails to yield a decrease in the value of the objective 
function, the pattern step is reti'aced and the search for pattern 
direction is made with half the step siae. This is carried on 
till the convergence is reached or the step siae becomes less than 
a specified value. 

3.4 APPLICATION TO ROLLING CC»4TACT 

As discussed in sec. 2.1, in rolling contact phenomenon, new 
nodes enter the contact aone and some leave it continuously, as 
rolling proceeds. It may be noted that the number of nodes in 
contact at any instant during steady rolling, remains constant. 
New nodes enter the contact aone at the leading edge and as the 
body keeps rolling, the locations of the nodes shift rearwards. 
Before the nodes leave contact, they have traversed through the 
complete contact aone. Tills is depicted in fig. 3.3, where nodes 
occupy different locations at time t^^, t^ and t^. During the 
rearward journey, the nodes undergo continuously varying 
stress-strain field along the contact length. Also, in rolling 
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motion, the spatial co-ordinates of every point of body changes, 
making modelling of such a problem difficult task. To overcome 
this, in present study, the concept of travelling finite elements 
is used. It jss may be noted that in case of steady rolling the 
stress-strain field and deformations over contact area remain same 
■with respecttsystem translating with speed of centre of roller. 

The method is based upon the fact that every material point 
or node on the periphery of rolling body passes through the entire 
contact field during its period of contact. For example, as shown 
in fig. 3.3, while the body rolls from instant t^ to node 1 
separates and node 5 enters the contact field. The external 
traction field is same for both situations. Under the steady 
state conditions, node 2 at instant t 2 , would have the same 
conditions of contact, as node 1 was having at the instant t . 
Similarly, node 3 replaces node 2 and so on. Tlius, the steady 
contact field is occupied by a newly formed set of nodes as the 
body rolls. 


In the present method, the above discussed feature of rolling 
contact is exploited. So, instead of actually rolling the body, 
the contact field conditions are shifted rearward over the contact 
nodal points. The contact zone is considered to be formed by the 
same set of nodal points for all iterations. 

The methodology of implementing the above algorithm is 
discussed herewith. To start with, the transient response at 
instant t is assumed to be known and rolling is assumed to take 
place in steps of one-degree which coincides with the distance 
between nodes in contact. 
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Let the contact be at n nodes as shown in fig. 3. 4 and the 
direction of rolling as shown. From the previous discussions the 
displacement boundary condition in tangential direction at the i*^^ 
node at instant t+At is taken as : 

u. , at time t+At = u...,, at time t 
1 1+1 

subscript i represents a node in contact zone. 

The normal displacement boundary conditions which are the 
results of kinematic constraints are taken constant for all 
iterations. Using these displacement boundary conditions, the 
response at instant t+At is obtained in terms of contact nodal 
forces. Now, the friction condition is checked at all the contact 
nodes to decide the contact status. Those nodes for which the 
ratio of tangential to normal traction exceeds the coefficient of 
friction, are allowed to slip. The extent of slip is determined 
using the miniDium dissipation of energy approach. The direction 
of opposes the slip at the node and the slip is equal to the 
change in displacement of that node while body has rolled from 
instant t to t+At. The equivalent of the slip in present method 
can be taken as ( ■ The subsequent iterations are 
carried out in similar manner till the response converges to a 
steady state one. 


The solution of rolling contact is obtained in two stages. 
In first stage, for the given external traction conditions, the 
static equilibrium solution is obtained. Then, in the second 
stage, the rolling criterion, as discussed above, is used to get 
the final steady state solution. The initial solution at any 
stage is obtained by taking sticking conditions at all nodes in 
contact. To decide about the region of slip, the nodes are 
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checked for the relative magnitude of taxigential and normal 
traction. At the nodes, satisfying the condition of slip, the 
extent of slip is determined using the principle of minimum 
dissipation of energy. The nodes are allowed to slip sequentially 
and their order of sequence is procedurally determined such that 
the resultant configuration has followed the path of minimvuB 
dissipation of energy. The friction condition is put at the nodes 
which are set to slip and the tangential displacements are 
obtained. 

3- 5 IMPLEMENT ATI <»f DETAILS 

In addition to putting the boundary conditions the 
solution of any problem requires the considerations of many 
other factors in order to make the program efficient and 
computationally economical. 

The global stiffness matrix [K3 , here, is symmetric and 
banded. To exploit these characteristics, it is stored in skyline 
form. Routines described in reference [Fellipa, 1974] are used 
for solution of matrices of such data structures. 


In this problem, the contact conditions i.e., a node coming 
into contact or a node separating or change in friction status of 
slipping or sticking when imposed, will change the stiffness 
matrix. Hence, the stiffness matrix gets modified in every 
iteration. However, the changes are only in the equations 
corresponding to the possible contact nodes. Also, the known 
displacements of the nodes can be eliminated from the solution 
system. Hence, a multi level siib-structuring scheme is made use 
so that a smaller matrix is handled in the iterative procedure. 
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The sub- structuring scheme Is necessary not only to save the 
computational time but also because the contact displacement and 
force relations, if imposed, render the stiffness matrix 
unsymmetrical , and hence requiring full storage. Sub-structuring 
enables smaller matrix, consisting of equations corresponding to 
the contact nodes only, to be stored and operated upon in every 
iteration. The main steps of sub- structuring are given in the 
Appendix A . 

To obtain the response for an incremental step, the chaxiged 
contact conditions are imposed iteratively, till the solution 
becomes consistent. The contact conditions, as mentioned in 
section 3.2, are imposed on the condensed stiffness matrix, which 
economizes the computer efforts. The exact methodology of 
imposition of these conditions is discussed in APPENDIX B. 
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At time t*) 


^2 


.t3 
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Fig.3.3 Different instances of rolling motion. 



1 2 i-1 i i+1 n 


Fig. 3.4 Typical representation of n ' nodes being 
in contact . 


CHAPTER IV 


RESULTS AlO DISCUSSION 


As discussed in earlier chapters two probleras are analysed to 
investigate the effectiveness of the pi'oposed algorithm. The 
first one is a plate with a hole having an oversized rigid 
inclusion, subjected to tjuasi- static loading and unloading. In 
second problem, the steady rolling of a cylinder on a rigid 
surface is studied. The validity of the results is also 
discussed. 


4.1. A Plate with a hole and inclusion 


As explained earlier the plate is assumed to be thick (Plane 
Strain condition) with a circular hole. An oversized rigid 
circular inclusion is embedded in the hole. The plate is 
subjected to a uniform farfield external loading (fig. 4.1a). The 
loading is assuD^ed to be slowly varying, cyclic in nature, 
satisfying the quasi-static assumption. A finite plate with a/r = 
6 is taken to simulate a large plate. Because of syiomietry it is 
enough to analyse one-quarter of the plate and the geometry is 
shown in fig. 4.1b along with finite element mesh. As shown in 
the above figure, there is a fine mesh near the hole boundary and 
ooai'ser mesh away. The actual mesh used for analysis is still 
finer with 360 elements and 207 nodes. 


In what follows the oversisedness of the rigid inclusion is 
expressed in terms of an equivalent initial precompression stress, 


T 

pre 


as 


T 

pre 


$ 


S E 
r (1 + n) 


C 4.1 ) 
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S. aw'iO'v-S ‘■Vr pv<s.<iov'>^.;pvac;,sion 

r radius of hole 

E modulus of elasticity 

Poisson's ratio 

The material properties and geometry of the plate are tahen 
as 

E - 30.0E06 Pa 
^ = 0.25 
r = 1 cm 
a - 6 cm 
o ~ 1 . 5E-04 cm 

The external loading in terms of displacemei^its is applied. 
It is assumed that there are no locked displacements initially and 
the tensile load is applied gradually. 

4.1,1 Comparison of tlie Solution Methods : with and witliout the 
energy dissipation term 

The quasi-static problem is analy.sed by two methods to 
evaluate the effectiveness of the dissipation of energy term. 
Thus, while determining the slip, no-slip and separation sones by 
ndnimising tlie objective functioii the solutions are obtained by 
two methods : 

I objective function, consists the error due to 

constraints, e only { This method is referred as 'without 
dissipation of energy’ ) 

II objective function consists the error due to constraints, 

« as well as the energy dissipation, ( This method 

is referred as ’with dissipation of energy’ ) 

Fig. 4.2 presents the contact sones as function of load T 
(non dimensionalized as while the load is applied in 

tension for a coefficient of friction, p,'oi 0.1 by the method I 
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(’without dissipation energy") and fig. 4.3 presents the results 
of the same proDleiu by the method II {’y?ith dissipation energy 
term’ ) . From the figures it is clear that with the increase in 
load, slip zone appears first and expajids. Later the nodes near 
90 begin to separate aiid spread inward. It can be observed that 
there is no difference in the' results by the two methods, though, 
first method requires interpolation of load for any node going 
from 3io-slip to slip state. Figs. 4.4 and 4.5 presents similar 
results for a coefficient of friction, p, of 0.4. Qualitatively, 
the behaviour is same as those for p=0. 1 except that no-slip zones 
are bigger. Because of the coarseness of the mesh the data are 
not available very much near to 2^ and 92^ , the present curves are 
extrapolated (shown by da.shed line) for tlie sake of completeness 
and understanding. 

Figs. 4.6 and 4.7 present the variation of contact zones with 
load wlien the load is decreased from the maximuai tensile load 
(referred to as ’unloading tension'). Results have been obtained 
for a coefficient of friction. by the method I and the- 
method II respectively. During this unloading process, the nodes 
tend to slip in the direction opposite to that while loading and 
is opposed by frictional forces. As unloading proceeds, for a 
fraction of load step the slip zone becoines no-slip. The nodes in 
no-slip zone go to slip gradually and the size of separation zone 
decreases and it changes into slip zone. It can be seen that the 
same behaviour is obtained by both the methods, Wliereas at the 
higher coefficient of friction, p=0.4 (figs. 4,8 and 4,9), the 
metiiods differ. It can be observed that in case of method I (fig. 
4.8) whei^ any node enters the contact zone from separation, it 
slips. Later on, the nodes near 90*^ go to no-slip zone. On the 
other hand, "method II (fig. 4.9) shows that the nodes near 90^ 
directly go to no-slip zone. After investigating the two methods, 
it was found that in case of method I, the nodes near 90* simply 
oscillates between slip axid ixo-slip conditions. It requires a 
finer nxesh to get correct solution with method I, whereas method 
II can tackle such situations effectively even with present mesh. 
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After completely anloading, it is loaded in compi'ession to 
the same maximum value as that in tension. The figs. 4.6 to 4.9 
show the behaviour during that process and found to be comparable 
for both methods. 

During unloading from maximum compression, in case of method 

I (’without dissipation of energy’), fig. 4.10 (for n = 0.1) shows 
that there is complete no-slip zone at the interface initially and 
slip zone appears as unloading continues. Wliereas, by the method 

II (’with dissipation of energy’) some of the nodes slip in the 
same direction as that while loading. As the unloading proceeds 
some of those Tiodes stick for some time and then start slipping iXi 
tiie reverse direction. This complex behaviour is not obtained by 
the method I procedure. Praobably, it may be possible to obtain 
this complex bahaviour by the method I also, if a very fine mesh 
and small load steps are used. Thus, it can be inferred that 
minimization of dissipation of energy during loading steps aids in 
the determination of the complex contact conditions more 
accurately with ease. 

4.1.2 Stress Variation at the Interface 

The distribution of normal and tangential stresses at the 
interface is presented in figs, 4.14 to 4.17, with stresses being 
noii-dimensionalised with respect to farfield applied stress, lor 
different coefficients of friction. 


Fig. 4.14 presents the distribution of normal stress at the 
maximum tensile load. As expected the stress is maximum near 0"^ . 
It can also be observed that, maximuBi, decreases with increase 
in coefficient of friction. There is a slight deviation near the 
plane of symmetry, which may be due to angular contact area 
[Gaertner, 1977] and use of lower order element in the region of 
high stress gradient. There exists a separation zone also at this 
load, shown as zero contact stresses. 
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The distribution of tangential stress at maximum tensile load 
is showii in fig. 4.15. The curves are extrapolated, shown by 
dashed lines, to represent values of stress near to and 

outermost edge of contact. The disti’ibution is consistent with 
friction conditions present at interface. 

Figs. 4.16 and 4.17 present normal and tangential stress 

distributioxi at the maximum compressive load. It may be seen that 

as the effect of friction increases the o becomes flatxer and the 

n 

maximum has a higher value for lower coefficients of friction 
wnereas mauiimuift is lower for lower coefficients of friction. 
The normal stress distribution has depression near both the plane 
of syiraiietry. Tlie effect of friction upoxi it i.s as expected. 


4.2 A Cylinder roiling on a Rigid plane suri acre 

As mentioned earlier in the case of a cylixider rolling ou a 
horizontal groujad, the rolling is assruned to take place at a 
uiiiiorm speed and speed is vt;ry low. As tlie horizontal ground 
IS assumed to be rigid aiid the cylixider is elastic, for xhe fiiiiXe 
element axialysis it is enough to discrete the cylinder. Details 
of the geometry and the typical mesh used is shown ixi fig. (4.18 a 
aiid b). Mesh is fixier xiear the contact zone with the ground. 
Actual mesh used is still finer with 1548 element-s and 847 nodes. 
Tiie nodes within fine mesh are positioned at one degree of axigle 
apart on the pei'iphery. The flat bed is modelled as an infiiiite 
dimensional bed but all the displacements of the bed are taken 
aero to depict, rigidity. Tlie W axid P_are normal arid driving 
taxigeiitial load on the cylixider and is considered for losses 
like windage, etc and in the present study it is taken constant as 

Fg/W = 0.02. 

As explained in sect. 3.4, the initial configuration is taken 
to be the static solution and rolling is effected in steps of one 
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degree so that the nodal contact conditions are translated to the 
left by one node at a time. In each iteration nodal conditions 
are modified to satisfy the contact conditions in a sequence so as 
to minimize the dissipation of energy. Tliis procedure is 
continued until steady contact conditions are obtaiiied arriving at 
steady rolling configuration. 

The solution obtained, satisfies all the physical laws 
applicable and all the characteristics of the rolling contact. 
Figures 4.19 and 4,20 show the distribution of normal and 
tangential forces at the contact nodes for different coefficients 
of friction, for a P/h’ = 0.04. The nodal force components are 
non-dimensionalised. The location of centre of tht cylinder is 
represented by a single dash centre line in all the figures. As 
discussed in various literatures [Moore, 1975], the slip regions 
exist over the leading and trailing edge only. The intez'mediate 
portion of the contact is of no-slip nature. Tlie exi-ent of slip 
over the leading edge{L.E. ) is very small(~ 1*^) and of significant 
magnitude over the trailing edgeCT.E. ; (~ 4-5°), thus agreeing 
with the concept of micro and macro slip over leading and trailing 
edge . 


The overall contact area and normal stress distribution are 
essentially same for different coefficients of friction (fig. 
4.19). As can be seen in the fig. 4.20, the tangential stress 
varies linearly near the middle portion of the contact and attains 
maximum v&lues near the edge arid assumes zero values at the ends 
of contact zone. The overall maximum occurs in the trailing edge 
and higher for higher coefficients of friction. Fig 4.21 shows 
the size of various contact zones at different coefficients of 
friction. It is clear from this figure that there occurs slip 
zones at the either edge ( As already mentioned ) and they 
decrease significantly as the coefficient of friction increases. 
Figure 4.22 and 4.23 presents the variation of normal and 
tangential nodal force components in the contact zones for various 
P/W values. This brings out the fact, as expected, that as the P/W 
increases the contact area shifts towards the leading edge. Tlie 



j-es£. 


remains same whereas 


normal st/ress disnric>u.t.ion more or 
nangent-iai increases with F,/n. 


hhe 


Effect of changing ? and W keeping P/W constant on the 
contact stresses can be »een i.n fig. 4.24 and 4.25. It can be 
observed that the overall contact length increases with P, but the 
variation is of the same nature. 


Energy dissipation per degree rotation is evaluated for 
different oases «nd presented in fig. 4.26 & 4.27. in fig. 4.26 
tne variation of energ’/ dissipation with P/W for higher 
coefficients friction is shown. It may be noted thcAt tlie 

eaergi dissipation increases with F/k’. Fig. 4.27 presents 
variation of ciiergs’ dissipatiOia with coefficient of friction for 


various. P/’W values. The dissipation of energy is low 
low 'I'uef f icient of friction and increases v/ith tfie n, 
mar.iiuuTi. value and decreases with further increase of w. 


at Cl very 
attains a 
Thus, it 


clearly shows tiiat for lowei dissipation of energy, 


it 1 .? required 


tC' have a very low or high value of coefficient of friction. 



Fig. 4.1 (a) A plate with a oversized rigid inclusion 



Fig. 4.1 (b) Representative mesh and prescribed boundory 
conditions . 
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Fig.4.2 Interface zones white loading in tension from 
undeformed state (without minimizing energy 
of dissipation) 



NODE POSITION, 9 (deg) 


3 



F'lQ A 3 Intcrfoc© zones while loodinQ intension from 
undeformed state (using minimum dissipation 
of energy approach ) 
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Fig.4.4 Interface zones while loading intension from 
undeformed state (without minimizing energy 
of dissipation) 
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Fig. 4.5 Interface zones while loading in tension from 
undeformed state, for )i = 0.4 (using 
minimum dissipation of energy approach) 



T /Tpre 


Fig.4.6 Interface zones while unloading tension. 

(without minimizing energy of dissipation) 
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Fig. 4.7 Interface zones while unloading tension and 
loading in compression (using minimum 
dissipation of energy approach) 



48 



Fig.4.8 Interface zones while unloading tension and 
loading in compression (without minimizing 
energy of dissipation) 
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Fig. 4.9 Interface zones while tension unloading and 
loading in compression , for p=0.4 (using 
minimum dissipation of energy approach) 




NODE POSITION, e (deg.) 
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Fig.4.10 Interface zones while unloading compre^Jon 
and loading in tension (without minimizing 
energy of dissipation) 
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Fig.A.11 Interface zones while unloading compression 
and loading in tension for p.=0.1 (using mini- 
mum dissipation of energy approach) 


] j30V;i 
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F‘ig.4.12 Interface zones while unloading compression 
and loading in tension (without minimizing 
energy of dissipation) 




NODE POSITION, 0 (deg.) 
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Fig. 4.13 Interface zones while compression unloading 
and loading intension, p. = 0.4 (using minimum 
dissipation of energy approach) 




NODE POSITION, 0 (deg.) 

Fia 4 14 Variation o1 normal stress over contact length 
■ for different n (at maximum tensile load, 

T/Tpre = 0.7) 
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NODE POSITION, 9 (deg.) 


Fig.4.16 Variation of normal contact stress over contact 
length, at different ji (at maximum compressive 
load, T/Tpre =-0.7) 







Fig.4.18(Q) Cylinder rolling on rigid surface. 



Fig.4.18(b) Representative mesh 
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Effect of friction upon normal force distribution. 
(W/ED r 0.834) 





COEFFICIENT OF FRICTION 


:l 



Fig.4.21 Variation of slip and no-slip zones over contact 
area with friction (at different P/W) 

(W/ED = 0.834) 
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T.E. NODES IN CONTACT L.E. 

Fig.4.22 Distribution of normal force over contact length, 
at different P (at constant p and W ) 

(W/ED =0.834) 



/ED (xIO 
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NODES IN CONTACT 


FigA23 Distribution of tangential force over contact 
length, at different P (at constant jj and W) 
(W/ED=0.834) 
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Fig.4.24 Distribution of normal contact force at constt 
P/W a u (W/ED =0.834) 
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Fig.4.25 Distribution of tangential contact force at 
constant P/W & H (W/ED =0.834) 






0 0.2 0 . 

Coefficient of Friction, ji 

Fig.4.27 Effect of friction upon energy of dissipation 
(Ed) at different P/W (W/ED = 0.834) 



CHAPTER V 


CONCLUSIONS AND SUGGESTIONS FOR FUTURE WORK 


The solutions of dynamic contact problems have been presented 
using the finite element and minimum dissipation of energy 
principle. Based upon the results presented in chapter- IV the 
following conclusions can be drawn : 

1- Using suggested procedure it is not required to keep 
track of the direction of frictional resistance during 
the various stages of loading. The procedure itself takes 
care of all such situations. 

2. In complex contact conditions the present method leads to 
more exact results with ease. For the same accuracy by 
the method without dissipation energy will require very 
fine mesh and very small load steps. 

Suggestions for Future Work : 

1 . More accurate element can be used to make the prograu^ 
more efficient. 

A more generalised rolling contact problem can be 
attempted wherein spin, creepage, etc. are also 
considered. The feasibility of employing this algoritlim 
in solving the three dimensional rolling problems can be 
explored. 


2 . 
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APPENDIX A 

SUB-STRUCTURI HG 

The equilibrium equation 
[ K ] { U > = I F 1 
is re-arranged as showii 


■ EK,] 

ty 

[K3] • 




- IF,1 ' 


[K4I 

[Kj] 



► 3 *1 

0 

[K33T 


tKg] _ 


< J 


{F.I 

J 


where, 

are pescribed displacements degrees of freedom (dof) 
are displacement dof corresponding to free nodes 
Do are contact dof and prescribed force dof 

Since {U,} is knov?n, the equation reduces to 


■ [K^I 

[Kj] ■ 


™2>'' 

i- zz •* 

riF„n 

c* 

. ^ i 

► 

. [ K 3 I ^ 

[Kg] . 


. <D3>. 




.'*'31^- 


where 

<^21^ = -I '^2' 

Since {Fg) corresponds to the internal dof, so 

{F 2 } = {0} 

By solving the equations for 

[K,y3 {n 3 } = {F 3 } + {F^} 



xin 


where 

[K^D = [Kg] - [Kg]*^ 

{F^} = {Jj} - [K5]''tK^]''-{F2j} 

The evaluation of [K^] and {F^} can be done in preparation 
stage. In the evaluation of these matrices no explicit inversion 
is done. Back substitution with previous factorised matrices is 
done to compute products like ( [K^] [Kg] , . . - , [Kg] 

are not stored separately but are obtained procedurally from the 
skyline stored stiffness Diatrix. [K.^] is stored in full storage 
mode. All the conditions on displacements {1)3} applied on 

[K.^]. After solving for {U^}, {Ug} dof can also be obtained by 
substituting for {U^} in appropriate equation. 
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APPENDIX B 

IMPOSITIOM OF CONTACT BOUNDARY CONDITIONS 


The conditions of slip and no-slip over the contact region, as 
mentioned in section 3.2, are reproduced here for convenience. 

For stick or no-slip: 
u - S 

n 1 1 

For slip: 


subscript 1,2 corresponds to nodes. These conditions are put 
in the condensed stiffness matrix. Their imposition on the matrix 
is explained by an example given below as: 

Let consider a condensed stiffness matrix [K^l for the dof 
corresponding to two nodes, one sticking and other slipping. The 
equation is given as : 


whex'e 


£K^3 { U } = { F } 



{ U } is nodal displacement vector and { F } is nodal force 


vector . 
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Let., the first two dof corresponds to sticking node and 
remaining two to slipping one. The [K^3 is modified while 
imposing the above conditions. To apply the friction condition at 
the slipping node, the ^ times 4^^^ row is added to the 3^'^ row, 
making right hand side zero. The sign of m is chosen such that 
opposes u^. The displacement conditions are imposed as shown: 
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d +AJd 
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VI 


s 

J n3! 
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The above equation now can be solved for displacements . 



